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Abstract 

Improved algorithms for the solution of the time- 
dependent Euler equations are presented for unsteady 
aerodynamic analysis involving unstructured dynamic meshes. 
The improvements have been developed recently to the spatial 
and temporal discretizations used by unstructured grid flow 
solvers. The spatial discretization involves a flux-split 
approach which is naturally dissipative and captures shock 
waves sharply with at most one grid point within the shock 
structure. The temporal discretization involves an implicit 
time-integration scheme using a Gauss-Seidel relaxation 
procedure which is computationally efficient for either steady 
or unsteady flow problems. For example, very large time 
steps may be used for rapid convergence to steady state, and 
the step size for unsteady cases may be selected for temporal 
accuracy rather than for numerical stability. Steady and 
unsteady flow results are presented for the NACA 0012 airfoil 
to demonstrate applications of the new Euler solvers. The 
unsteady results were obtained for the airfoil pitching 
harmonically about the quarter chord. The resulting 
instantaneous pressure distributions and lift and moment 
coefficients during a cycle of motion compare well with 
experimental data. The paper presents a description of the 
Euler solvers along with results and comparisons which assess 
the capability. 


Introduction 

Considerable progress has been made over the past two 
decades on developing computational fluid dynamics (CFD) 
methods for aerodynamic analysis.1.2 Recent work in CFD has 
focused primarily on developing algorithms for the solution of 
the Euler and Navier-Stokes equations. For unsteady 
aerodynamic and aeroelastic analysis, these methods generally 
require that the mesh move to conform to the instantaneous 
position of the moving or deforming body under consideration. 
Many of the methods that are currently being developed assume 
that the mesh moves rigidly or that the mesh shears as the 
body deforms. These assumptions consequently limit the 
applicability of the procedures to rigid-body motions or 
small-amplitude deformations. Furthermore, these methods 
of solution typically assume that the computational grid has an 
underlying geometrical structure. As an alternative, 
algorithms have been developed recently which make use of 
unstructured grids. 3-12 | n j wo dimensions these grids are 
typically made up of triangles, and in three dimensions they 
consist of an assemblage of tetrahedra. The unstructured grid 
methods have distinct advantages over structured grid methods 
in that they can easily treat the most complex of geometric 
configurations as well as flow conditions, and that the 
unstructured grid can be moved to treat realistic motions and 
structural deformations of these configurations. 1 0-1 2 

The results presented by the author in Refs. 11 and 12 
demonstrated that (1) the methods produce solutions of 
comparable accuracy to results obtained using structured grid 
flow solvers. 11 and (2) that the unstructured grid 
methodology can easily analyze complex aircraft geometries 
undergoing structural deformation. 12 The methods of Refs. 11 
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and 12, however, use a spatial discretization based on central 
differencing with explicit artificial dissipation, and use a 
temporal discretization involving explicit time-marching 
based on a multi-stage Runge-Kutta time integration. The 
explicit artificial dissipation used in such schemes tends to 
smear shock waves over several grid cells and requires the 
tuning of free parameters that scale the dissipation. Also, the 
explicit Runge-Kutta time-integration has a step size that is 
limited by the Courant-Fredricks-Lewy (CFL) condition to 
very small values. Consequently, thousands (and occasionally 
tens of thousands) of time steps are required to obtain steady- 
state solutions, and thousands of steps per cycle of motion are 
required for unsteady solutions. Therefore, the purpose of the 
paper is to report on improvements that have been developed 
recently to the spatial and temporal discretizations of the 
unstructured grid flow solvers which resolve the numerical 
issues described above. The spatial discretization now 
involves a so-called flux-split approach, which is similar to 
discretizations presented in Refs. 10. 13. and 14 based on 
either the flux-vector splitting (FVS) of van Leer 1 5 or the 
flux-difference splitting (FDS) of Roe. 1 6 These flux-split 
discretizations account for the local wave-propagation 
characteristics of the flow and they capture shock waves 
sharply with at most one grid point within the shock 
structure. A further advantage is that these discretizations 
are naturally dissipative and consequently do not require 
additional artificial dissipation terms or the adjustment of 
free parameters to control the dissipation. Furthermore, the 
temporal discretization has been changed to an implicit time- 
integration scheme involving a Gauss-Seidel relaxation 
procedure similar to discretizations presented in Refs. 17 and 
18. This relaxation scheme is unconditionally stable and thus 
allows the selection of the step size based on the temporal 
accuracy dictated by the problem being considered, rather than 
on the numerical stability of the algorithm. Consequently, 
very large time steps may be used for rapid convergence to 
steady state, and an appropriate step size may be selected for 
unsteady cases, independent of numerical stability issues. 
Steady and unsteady results are presented for the NACA 0012 
airfoil to demonstrate applications of the new Euler solvers. 
The unsteady flow results were obtained for the airfoil 
pitching harmonically about the quarter chord. The paper 
presents a description of the Euler solvers along with results 
and comparisons which assess the capability. 


Euler Equations 

In the present study, the flow is assumed to be governed by 
the two-dimensional time-dependent Euler equations which 
may be written in integral form as 

in J/Qdxdy + J (Fdy-Gdx) = 0 m 

a an 

where the vector of conserved variables Q and the convective 
fluxes F and G are given by 
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The contravariant velocities U and V are defined by 

U=u-x t V=v-y t ( 3 ) 

where x t and y t are the grid speeds in the x and y directions, 
respectively, and the pressure p is given by the equation of 
state for a perfect gas 


P = (Y -1) t e - ^ P ( u 2 +v 2 ) ] (4) 

The above equations have been nondimensionalized by the 
freestream density and the freestream speed of sound a*,. 


and U and V are the corresponding contravariant velocities 
u„ ( „-x, <e.) 
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The flux vector H is split in a one-dimensional fashion into 
forward (H + ) and backward (H~) vectors for |U| £ a as 


where 
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Flux-Vector Splitting 


The spatial discretization based on flux-vector splitting is 
a cell-centered scheme where the flow variables are stored at 
the centroid of each triangle and the control volume is simply 
the triangle itself. The boundary integral of Eq. (1) is 
approximated by using the flux-vector splitting of van 
Leer , 1 8 | n this method the flux vectors are split into forward 
and backward contributions which are continuously 
differentiable even at sonic and stagnation points. The scheme 
is derived as follows. For each edge of a given triangle, the 
fluxes are first rotated into a locally Cartesian coordinate 
system x-y with the principal direction being perpendicular 
to the edge. The flux in this direction is defined as 


HAs = T(FAy - GAx) = < 
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pUu + p | 
pUv 

eU + puj 
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where the transformation matrix T is given by 
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In Eqs. (5) and (6), Ax and Ay are the directed lengths of the 
edge in the x and y coordinate directions, respectively, and 

As 2 = Ax 2 + Ay 2 . Also, u and v are the Cartesian velocity 
components perpendicular and parallel to the edge defined by 
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The resulting split fluxes are finally rotated back, into the 
original coordinate system so that 

FAy-GAx = T-'[H + (q-) + H-(q*)] (12) 

where the notation H+(q*) and H*(q+) indicates that the 

fluxes H are evaluated using upwind-biased interpolations of 
the primitive variables q. For a given triangle j, for example, 
and considering the diagram in Fig. 1(a), the upwind-biased 
interpolation for q' along the edge between triangles j and k is 
defined by 

q _ = q, +7K1-K) a.+(i+k) aj 

where 

A + = q k -qj 
A- = q, -q, 

In Eqs. (13) and (14), q^ and q k are the vectors of 
primitive variables at the centroids of triangles j and k, 
respectively, and q jf the vector of primitive variables at 

node i, is determined by an average of the flow variables in the 
triangles surrounding node i. The upwind-biased interpolation 
for q + along this edge is determined similarly using the flow 
variables at centroids j and k and the flow variables at node I. 
The parameter k in Eq. (13) controls a family of difference 
schemes by appropriately weighting a. and A^. On 
structured meshes it is easy to show that k - - 1 yields a fully 
upwind scheme, k « 0 yields Fromm's scheme, and k - 1 yields 
central differencing. The value k - 1/3 leads to a third- 
order-accurate upwind-biased scheme, although third-order 
accuracy is strictly correct only for one-dimensional 
calculations. Nevertheless, k - 1/3 was used in the 
calculations presented herein. 


(13) 

(14a) 

(14b) 
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On highly stretched meshes, the formula for A + | s 
modified to be 


A = 


TTb (q K-V 


(15) 


<r - Qj + J [ (1 -KS)A_+ (1 + KS)A + ] (17) 

where s Is the flux limiter. In the present study, the 
continuously differentiable flux limiter of Ref. 21 was 
employed which is defined by 


where a and b are the distances from the midpoint of the edge to s = + 

the centroids of triangles j and k, respectively, as shown in A 2 + a 2 +e (18) 

Fig. 1(b). This formula weights the flow variables in the . . (l “ + 

interpolation formula (Eq. (13)) differently to account for wnere c a very small number to prevent division by zero in 

the stretching of the mesh. For example, by substituting Eq. smooth regions of the flow. 

(15) into Eq. (13) and letting k - 1 yields 


q" = — - — q + — — q 
a+b M i a + b H k 


(16) 


For the case shown in Fig. 1(b), Eq. (16) clearly gives more 
weight in the calculation of q* to the flow variables at centroid 
j than to the flow variables at centroid k, since b>a. 
Furthermore, in calculations involving upwind-biased 
schemes, oscillations in the solution near shock waves are 
expected to occur. To eliminate these oscillations flux limiting 
is usually required. The flux limiter modifies the upwind- 
biased interpolations for q* and q + such that, for example 


• centroid 



( a ) centroids and nodes used in construction 
of upwind-biased flow variables. 


Flux-Difference Splitting 

The spatial discretization based on flux-difference 
splitting is a cell-centered scheme where the flow variables 
are stored at the centroid of each triangle and the control 
volume is simply the triangle itself. The flux balance becomes 

XT-'HAs = *£T- l [H(q-) + H(q + )]As 

-IX T-'fAKcr-cr) as d9) 


where the first term on the right-hand side is simply the 
average of the original flux H evaluated using the flow 
variables on each side of the edge. The second term on the 
right-hand side represents a flux difference since it involves 
the product of the difference in flow variables (Q+- Q’) and 
the flux jacobian A which is defined as 



In this context the flux jacobian can be rewritten by a 
similarity transformation as 

A = RAR _1 (21) 

where A is a diagonal matrix whose diagonal elements are the 
characteristic speeds U, U, U+a, and U-a. The variable U 
is the contravariant velocity normal to the edge being 

considered and a is the local speed of sound. The notation IAI 
indicates that the flux jacobian is evaluated by taking the 
absolute value of the characteristic speeds and by using so- 
called Roe-averaged16 flow variables (indicated by the tilde). 



( b ) distances between centroids and midpoint 
of edge used in Eqs. (15) and (16). 

Fig. 1 Diagrams illustrating details of flux-split Euler 
algorithm implementation. 


In the FDS scheme, the conserved flow variables on either 
side of a given edge, Q+and Q-, are determined by first 
calculating the upwind-biased primitive flow variables q’ and 
q+ and then converting them from primitive to conserved. The 
interpolation formula for q*. for example, is identical to Eq. 
(17) including the flux limiter defined by Eq. (18). 
Therefore the same parameter k in Eq. (17) controls a family 
of difference schemes, ranging from fully upwind to central 
differencing, by appropriately weighting A _ and A + as in the 
FVS scheme. The calculations presented herein used k « 1/3, 

implicit Temporal Discroti?atinn 

The Implicit relaxation algorithm is formulated by first 
approximating the time derivative in the Euler equations by 

.30. 2+4. aQ 2+» Q‘-Q n » Q n -Q n ~ 1 

dt 2 At 2 At 2 At 22 ^ 

n +i * 

where AQ=Q -Q and where the parameter $ controls 
the temporal order of accuracy. For example, the scheme is 
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first-order-accurate in time if $-0 and the scheme is second- 
order-accurate in time if $*1. For an implicit temporal 
discretization, the flux H must be treated at time level (n+1) 
which is accomplished by linearizing according to 


H n+1 =H’ 

3Q 


dH 


Q = Q 


. AQ 


(23) 


where 3H/3Q is the flux jacobian A as discussed above. Also, 
Eqs. (22) and (23) involve Q*, the vector of flow variables at 
an iterate level O which is normally taken to be time level 
(n). For unsteady applications, however, subiterations may 

be performed to drive Q to Q n+1 and thus minimize 
linearization and relaxation errors. 


except of course, the residual is computed using FDS. The 
second way is to construct approximate jacobians using Eq. 
(21) and the fact that the forward and backward jacobians 
should have non-negative and non-positive eigenvalues 
(characteristic speeds), respectively. This is done to produce 
a diagonally dominant system of equations for numerical 
stability and is accomplished by defining 


a + =ra" r' 1 

where 

A ♦ A + I A | 

A — — 


a' = ra-r _1 


- A -|A| 
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(26) 


(27) 


For FVS the forward and backward fluxes are linearized 
for a given triangle j as 

(q“) + H” (q*) ]“*' As 

= £t-'[H + (q-)+H-(q*)]*As 

+ ( X T~' A+ As ] AQ; + £ T-'A-AsACL 

(24) 

In this equation the last summation on the right hand side 
involves A Q m , the change in the flow variables in the three 
triangles adjacent to triangle j. Also, the exact jacobians 
A + and A‘ are determined by differentiation of H + and H‘ by 
the conserved variables Q. By combining Eqs. (22) and (24), 
the Euler equations are discretized as 


The results presented herein were obtained using the first of 
these two approaches. 


Results and Discussion 

To assess the new Euler solvers, calculations were 
performed for the NACA 0012 airfoil. These results were 
obtained using the unstructured grid shown in Fig. 2 which 
was generated using the advancing front method. 6 * 19 The grid 
has 3300 nodes, 6466 triangles, and extends 20 chordlengths 
from the airfoil with a circular outer boundary. Also there 
are 110 points that lie on the airfoil surface. This is the same 
mesh that was used to obtain the results that were presented in 
Ref. 11. Steady-state calculations were performed for the 

airfoil at a freestream Mach number of M_ = 0.8 and an angle 
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where I is the identity matrix and "area" is the area of 
triangle j. Direct solution of the system of simultaneous 
equations which results from application of Eq. (25) for all 
triangles in the mesh, requires the inversion of a large matrix 
with large bandwidth which is computationally expensive. 
Instead, a Gauss-Seidel relaxation approach is used to solve the 
equations whereby the summation involving A Q m is moved to 
the right hand side of Eq. (25). The terms in this summation 
are then evaluated for a given time step using the most 
recently computed values for A Q m . The solution procedure 
then involves only the inversion of a 4 x 4 matrix 
(represented by the terms in square brackets on the left hand 
side of Eq. (25)) for each triangle in the mesh. Also, although 
the procedure is implemented for application on (randomly- 
ordered) unstructured meshes, it is not a point Gauss-Seidei 
procedure. The method is in fact more like line Gauss-Seidel 
since the list of triangles that make up the unstructured mesh 
is re-ordered from upstream to downstream, and the solution 
is obtained by sweeping two times through the mesh as dictated 
by stability considerations. The first sweep is performed in 
the direction from upstream to downstream and the second 
sweep is from downstream to upstream. For purely 
supersonic flows the second sweep is unnecessary. 


For FDS the exact jacobian A is too expensive to compute 
and thus an approximate jacobian is normally used. There are 
several ways to accomplish this, two of which are described as 
follows. The first way is to simply use the forward and 
backward jacobians from the FVS scheme as in Eq. (25), 


of attack of a 0 =1.25°. Unsteady calculations were performed 
for the airfoil pitching harmonically about the quarter chord 
with an amplitude of =2.51° and a reduced frequency based 

on semichord of k * 0.0814 at M^ = 0.755 and a Q = 0.016°. 

These calculations are compared with the experimental data of 
Ref. 20. 

Steady Flow Results 

Steady flow results were obtained for the NACA 0012 
airfoil using both the implicit time-marching of the present 
study and the explicit four-stage Runge-Kutta time-marching 
of Ref. 11. Results are presented first using flux-vector 



Fig. 2 Partial view of unstructured grid of triangles about 
the NACA 0012 airfoil. 
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splitting and then using flux-difference splitting for the 
spatial discretization. The explicit time-marching results 
were obtained using a CFL number of 2.5 (since the CFL limit 
is approximately 2.8) and the implicit time-marching results 
were obtained using a CFL number of 100,000. Such a large 
value was used for the implicit results since the relaxation 
scheme has maximum damping and hence fastest convergence 
for very large time steps. This is in contrast with implicit 
approximate factorization schemes which have maximum 
damping for CFL numbers on the order of 10. 

A comparison of the convergence histories between explicit 
and implicit time-marching for flux-vector splitting is shown 
in Fig. 3(a). The "error - in the solution was taken to be the 
L2 norm of the density residual. As shown in Fig. 3(a), the 
explicit solution is very slow to converge. This solution takes 
approximately 10,000 time steps to become converged to 
engineering accuracy, which is taken to be a four order of 
magnitude reduction in solution error. In contrast, the 
implicit solution is converged to four orders of magnitude in 
only approximately 500 steps and is converged to machine 
zero in less than 2000 steps. The implicit solution costs 
approximately 75% more per time step than the explicit 
solution because of the increased number of operations 
required to evaluate the flux jacobians. This increase in CPU 
time is far out-weighed by the faster convergence to steady 
state in that a converged solution is obtained with the implicit 
relaxation scheme with an order of magnitude less CPU time 
than the explicit scheme. The resulting steady pressure 


distribution is shown in Fig. 3(b). For this case there is a 
relatively strong shock wave on the upper surface of the 
airfoil near 62% chord and a relatively weak shock wave on 
the lower surface near 30% chord. The pressure 
distributions Indicate that there is only one grid point within 
the shock structure, on either the upper or lower surface of 
the airfoil, due to the sharp shock capturing ability of flux- 
vector splitting. Furthermore, the steady pressure results of 
Fig. 3(b) are of comparable accuracy in comparison with the 
numerous published results for this case such as those 
reported in Ref. 21. 

A comparison of the convergence histories between explicit 
and implicit time marching for flux-difference splitting is 
shown in Fig. 4(a). Similar to the solutions obtained using 
flux-vector splitting, the explicit solution here is very slow 
to converge. However, the implicit solution is again converged 
to four orders of magitude in only approximately 500 steps 
and is converged to machine zero is less than 2000 steps. 
These solutions, with either implicit or explicit time- 
marching, cost approximately the same as the corresponding 
solutions involving flux-vector splitting. The resulting steady 
pressure distribution is shown in Fig. 4(b). The pressure 
distribution again indicates that there is only one grid point 
within the shock structures due to the flux-difference 
splitting, and the shocks appear to be slightly more sharply 
captured in comparison with the shocks from the solution 
obtained using flux-vector splitting. This is because the FDS 
scheme has less dissipation than the FVS scheme. 
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(a) convergence histories. 


(a) convergence histories. 



(b) steady pressure distribution. 

Fig. 3 Comparison of steady-state results for the NACA 
0012 airfoil at M--0.8 and a 0 - 1-25° computed 
using flux-vector splitting. 



x/c 

(b) steady pressure distribution. 

Fig. 4 Comparison of steady-state results for the NACA 
0012 airfoil at M««0.8 and a 0 - 1.25° computed 
using flux-difference splitting. 
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Unsteady F low Results 

Unsteady results were obtained for the pitching NACA 
0012 airfoil using 250, 1000, and 2500 steps per cycle of 
motion with the Implicit time-marching and flux-vector 
splitting to determine the appropriate step size to ensure 
temporal accuracy for this case. Three cycles of motion were 
computed to obtain periodic solutions. The effects of step size 
on the instantaneous pressure distribution at kx-69° in the 
third cycle, which corresponds to an instantaneous pitch angle 
of a(t) » 2.34°, are shown in Fig. 5. This angle was selected 
for this assessment since it lies in the most sensitive part of 
the cycle. The results of Fig. 5 clearly indicate that with as 
few as 250 steps per cycle the upper and lower surface shocks 
both have inaccurate strength and location in comparison with 
the results obtained using 2500 steps per cycle. The results 
of Fig. 5 also indicate that the appropriate step size for this 
case is a time step corresponding to between 1000 and 2500 
steps per cycle of motion. This finding is consistent with the 
temporal refinement study of Ref. 22, where similar results 
were obtained using implicit approximate factorization 
solutions of the transonic small-disturbance and Euler 
equations. It is noted, however, that for easier cases involving 
higher reduced frequencies and smaller amplitudes of motion, 
as few as two or three hundred steps per cycle of motion are 
sufficient for temporal accuracy. 


The effects of performing subiterations per time step on 
the instantaneous pressure distributions at k x - 69° in the 
cycle corresponding to an instantaneous pitch angle of a(t) - 
2.34° are shown in Fig. 6. The calculations were performed 
using 250 steps per cycle of motion with the implicit time- 
marching and flux-vector splitting and parallel results were 
obtained using 0, 5, and 10 subiterations per time step. As 
discussed previously, the purpose of the subiterations is to 
minimize linearization and relaxation errors, similar to that 
which is done with approximate factorization schemes to 
minimize linearization and factorization errors. 22 As shown 
in Fig. 6, as the number of subiterations is increased the 
inaccuracies in shock strength and location are decreased. 
With 10 subiterations, for example, the instantaneous 
pressure distributions resemble closely those of Fig. 5, 
obtained using 2500 steps per cycle of motion and no 
subiterations. There is therefore a compromise between 
running large time steps with subiterations and running 
smaller time steps with no subiterations, since the CPU time 
is approximately the same. 

Instantaneous pressure distributions at eight points in 
time during the third cycle of motion from the 2500 steps per 
cycle solution using flux-vector splitting are shown in Fig. 7 
for comparison with the experimental data. In each pressure 
plot the instantaneous pitch angle a(x) and the angular 


cc(t) = 2.34° 
kx - 69° 



Fig. 5 Effects of step size on the instantaneous pressure distribution at kx * 69° and a(x)- 

2.34° during the third cycle of motion for the NACA 0012 airfoil pitching at 
» 0.755, a 0 - 0.016°, ai ■ 2.51°, and k = 0.0814 computed using flux-vector 
splitting. 


o(t) « 2.34° 
kx - 69° 



250 steps/cycle 



Fig. 6 Effects of performing subiterations per time step on the instantaneous pressure 

distribution at kx - 69° and a(x)» 2.34° during the third cycle of motion for the NACA 
0012 airfoil pitching at Moo « 0.755, cc 0 - 0.016°, ai * 2.51°, and k - 0.0814 
computed using flux-vector splitting. 
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Fig. 7 Comparison of instantaneous pressure distributions for the NACA 0012 airfoil 

pitching at - 0.755, a 0 • 0.01 6°, ai » 2.51 ", and k « 0.0814 computed using 
flux-vector splitting. 

Flux-difference splitting o Upper surface - Experiment 

□ Lower surface - Experiment 








position in the cycle kx are noted. During the first part of the 
cycle there is a shock wave on the upper surface of the airfoil, 
and the flow over the lower surface is predominately 
subcritical. During the latter part of the cycle the flow about 
the upper surface is subcritical, and a shock forms along the 
lower surface. The pressure distributions indicate that the 
shock position oscillates over approximately 25% of the chord 
along each surface, and in general, that the two sets of 
calculated results compare well with the experimental data. 
Similar to the steady flow results, the shock waves are 
captured sharply with at most one grid point within the shock 
structure. The calculated results, however, show the expected 
symmetry in the flow, in that the upper surface pressure 
distribution during the first half of the cycle is very similar 
to the lower surface pressure distribution during the second 
half of the cycle. The experimental data therefore appears to 
have been obtained at a slightly higher effective steady-state 
angle of attack than that reported in Ref. 20. Furthermore, 
the unstructured grid results of Fig. 7 are of comparable 
accuracy in comparison with published results obtained using 
structured grid methods for this case, such as those reported 
in Ref. 22. 

Similar comparisons between calculated and experimental 
instantaneous pressure distributions at the same eight points 
in time during the cycle are shown in Fig. 8. The calculated 


Euler 



a (deg.) 


(a) lift coefficient. 

Euler 



a (deg.) 


(b) moment coefficient. 

Fig. 9 Comparisons of coefficient versus instantaneous angle 
of attack for the NACA 0012 airfoil pitching at 
= 0.755, a 0 « 0.016°, cti » 2.51°. and k - 
0.0814 


results were obtained using the implicit time-marching with 
2500 steps per cycle of motion and the flux-difference 
splitting for the spatial discretization. The FDS pressures 
show similar features as the FVS pressures of Fig. 7 in that 
the shock waves are sharply captured within only one grid 
point within the shock structure. In general, the FDS 
pressure distributions also agree well with the experimental 
data. 

Comparisons of calculated and experimental lift and 
moment coefficients versus the instantaneous angle of attack 
are presented in Fig. 8. The lift coefficient is shown in Fig. 
8(a), and the moment coefficient is shown in Fig. 8(b). These 
coefficients show the variation as a function of angle of attack 
during a cycle of motion, and in general, the two sets of 
calculated results compare well with the experimental data. 
The comparisons of lift coefficient further indicate that the 
data was probably obtained at a higher effective steady-state 
angle of attack, since the experimental values are higher than 
the calculated values. Also, the largest difference between FDS 
and FVS coefficients which occur in the moment coefficient 
(Fig. 8(b)) are due to the sensitivity of the moment since the 
moment center is at the quarter-chord. The two calculated 
moment coefficients are not symmetric about one another 

because of the small angle of attack (a o - 0.016°) for this 
case. 

Concluding Remarks 

Improved algorithms for the solution of the time- 
dependent Euler equations were presented for unsteady 
aerodynamic analysis involving unstructured dynamic meshes. 
The improvements have been developed recently to the spatial 
and temporal discretizations used by unstructured grid flow 
solvers. The improved spatial discretization involves a flux- 
split approach which is naturally dissipative and captures 
shock waves sharply with at most one grid point within the 
shock structure. The improved temporal discretization 
involves an implicit time-integration scheme using a Gauss- 
Seidel relaxation procedure which is computationally efficient 
for either steady or unsteady flow problems. For example, 
very large time steps may be used for rapid convergence to 
steady state, and the step size for unsteady cases may be 
selected for temporal accuracy rather than for numerical 
stability. 

Steady and unsteady flow results were presented for the 
NACA 0012 airfoil to demonstrate applications of the new 
Euler solvers. The steady results showed that rapid 
convergence to steady state could be achieved with the implicit 
time-marching in comparison with results obtained using 
explicit time-marching. A factor of ten reduction in 
computational cost was obtained for the case that was 
presented. The unsteady results were obtained for the airfoil 
pitching harmonically about the quarter chord. The effects of 
step size and of performing subiterations per time step on the 
instantaneous pressure distributions during a cycle of pitching 
motion were demonstrated. These results indicated that the 
scheme was numerically stable for large time steps although 
smaller time steps were required to maintain temporal 
accuracy for the unsteady case that was considered. Also, the 
calculated instantaneous pressure distributions and lift and 
moment coefficients during a cycle of motion compared well 
with experimental data. 
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